home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / cher2k.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  13.0 KB  |  368 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE CHER2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        UPLO, TRANS
  8.       INTEGER            N, K, LDA, LDB, LDC
  9.       REAL               BETA
  10.       COMPLEX            ALPHA
  11. *     .. Array Arguments ..
  12.       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * )
  13. *     ..
  14. *
  15. *  Purpose
  16. *  =======
  17. *
  18. *  CHER2K  performs one of the hermitian rank 2k operations
  19. *
  20. *     C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C,
  21. *
  22. *  or
  23. *
  24. *     C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C,
  25. *
  26. *  where  alpha and beta  are scalars with  beta  real,  C is an  n by n
  27. *  hermitian matrix and  A and B  are  n by k matrices in the first case
  28. *  and  k by n  matrices in the second case.
  29. *
  30. *  Parameters
  31. *  ==========
  32. *
  33. *  UPLO   - CHARACTER*1.
  34. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  35. *           triangular  part  of the  array  C  is to be  referenced  as
  36. *           follows:
  37. *
  38. *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  39. *                                  is to be referenced.
  40. *
  41. *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  42. *                                  is to be referenced.
  43. *
  44. *           Unchanged on exit.
  45. *
  46. *  TRANS  - CHARACTER*1.
  47. *           On entry,  TRANS  specifies the operation to be performed as
  48. *           follows:
  49. *
  50. *              TRANS = 'N' or 'n'    C := alpha*A*conjg( B' )          +
  51. *                                         conjg( alpha )*B*conjg( A' ) +
  52. *                                         beta*C.
  53. *
  54. *              TRANS = 'C' or 'c'    C := alpha*conjg( A' )*B          +
  55. *                                         conjg( alpha )*conjg( B' )*A +
  56. *                                         beta*C.
  57. *
  58. *           Unchanged on exit.
  59. *
  60. *  N      - INTEGER.
  61. *           On entry,  N specifies the order of the matrix C.  N must be
  62. *           at least zero.
  63. *           Unchanged on exit.
  64. *
  65. *  K      - INTEGER.
  66. *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  67. *           of  columns  of the  matrices  A and B,  and on  entry  with
  68. *           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
  69. *           matrices  A and B.  K must be at least zero.
  70. *           Unchanged on exit.
  71. *
  72. *  ALPHA  - COMPLEX         .
  73. *           On entry, ALPHA specifies the scalar alpha.
  74. *           Unchanged on exit.
  75. *
  76. *  A      - COMPLEX          array of DIMENSION ( LDA, ka ), where ka is
  77. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  78. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  79. *           part of the array  A  must contain the matrix  A,  otherwise
  80. *           the leading  k by n  part of the array  A  must contain  the
  81. *           matrix A.
  82. *           Unchanged on exit.
  83. *
  84. *  LDA    - INTEGER.
  85. *           On entry, LDA specifies the first dimension of A as declared
  86. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  87. *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  88. *           be at least  max( 1, k ).
  89. *           Unchanged on exit.
  90. *
  91. *  B      - COMPLEX          array of DIMENSION ( LDB, kb ), where kb is
  92. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  93. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  94. *           part of the array  B  must contain the matrix  B,  otherwise
  95. *           the leading  k by n  part of the array  B  must contain  the
  96. *           matrix B.
  97. *           Unchanged on exit.
  98. *
  99. *  LDB    - INTEGER.
  100. *           On entry, LDB specifies the first dimension of B as declared
  101. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  102. *           then  LDB must be at least  max( 1, n ), otherwise  LDB must
  103. *           be at least  max( 1, k ).
  104. *           Unchanged on exit.
  105. *
  106. *  BETA   - REAL            .
  107. *           On entry, BETA specifies the scalar beta.
  108. *           Unchanged on exit.
  109. *
  110. *  C      - COMPLEX          array of DIMENSION ( LDC, n ).
  111. *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  112. *           upper triangular part of the array C must contain the upper
  113. *           triangular part  of the  hermitian matrix  and the strictly
  114. *           lower triangular part of C is not referenced.  On exit, the
  115. *           upper triangular part of the array  C is overwritten by the
  116. *           upper triangular part of the updated matrix.
  117. *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  118. *           lower triangular part of the array C must contain the lower
  119. *           triangular part  of the  hermitian matrix  and the strictly
  120. *           upper triangular part of C is not referenced.  On exit, the
  121. *           lower triangular part of the array  C is overwritten by the
  122. *           lower triangular part of the updated matrix.
  123. *           Note that the imaginary parts of the diagonal elements need
  124. *           not be set,  they are assumed to be zero,  and on exit they
  125. *           are set to zero.
  126. *
  127. *  LDC    - INTEGER.
  128. *           On entry, LDC specifies the first dimension of C as declared
  129. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  130. *           max( 1, n ).
  131. *           Unchanged on exit.
  132. *
  133. *
  134. *  Level 3 Blas routine.
  135. *
  136. *  -- Written on 8-February-1989.
  137. *     Jack Dongarra, Argonne National Laboratory.
  138. *     Iain Duff, AERE Harwell.
  139. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  140. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  141. *
  142. *
  143. *     .. External Functions ..
  144.       LOGICAL            LSAME
  145.       EXTERNAL           LSAME
  146. *     .. External Subroutines ..
  147.       EXTERNAL           XERBLA
  148. *     .. Intrinsic Functions ..
  149.       INTRINSIC          CONJG, MAX, REAL
  150. *     .. Local Scalars ..
  151.       LOGICAL            UPPER
  152.       INTEGER            I, INFO, J, L, NROWA
  153.       COMPLEX            TEMP1, TEMP2
  154. *     .. Parameters ..
  155.       REAL               ONE
  156.       PARAMETER        ( ONE  = 1.0E+0 )
  157.       COMPLEX            ZERO
  158.       PARAMETER        ( ZERO = ( 0.0E+0, 0.0E+0 ) )
  159. *     ..
  160. *     .. Executable Statements ..
  161. *
  162. *     Test the input parameters.
  163. *
  164.       IF( LSAME( TRANS, 'N' ) )THEN
  165.          NROWA = N
  166.       ELSE
  167.          NROWA = K
  168.       END IF
  169.       UPPER = LSAME( UPLO, 'U' )
  170. *
  171.       INFO = 0
  172.       IF(      ( .NOT.UPPER               ).AND.
  173.      $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
  174.          INFO = 1
  175.       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
  176.      $         ( .NOT.LSAME( TRANS, 'C' ) )      )THEN
  177.          INFO = 2
  178.       ELSE IF( N  .LT.0               )THEN
  179.          INFO = 3
  180.       ELSE IF( K  .LT.0               )THEN
  181.          INFO = 4
  182.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  183.          INFO = 7
  184.       ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN
  185.          INFO = 9
  186.       ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
  187.          INFO = 12
  188.       END IF
  189.       IF( INFO.NE.0 )THEN
  190.          CALL XERBLA( 'CHER2K', INFO )
  191.          RETURN
  192.       END IF
  193. *
  194. *     Quick return if possible.
  195. *
  196.       IF( ( N.EQ.0 ).OR.
  197.      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  198.      $   RETURN
  199. *
  200. *     And when  alpha.eq.zero.
  201. *
  202.       IF( ALPHA.EQ.ZERO )THEN
  203.          IF( UPPER )THEN
  204.             IF( BETA.EQ.REAL( ZERO ) )THEN
  205.                DO 20, J = 1, N
  206.                   DO 10, I = 1, J
  207.                      C( I, J ) = ZERO
  208.    10             CONTINUE
  209.    20          CONTINUE
  210.             ELSE
  211.                DO 40, J = 1, N
  212.                   DO 30, I = 1, J - 1
  213.                      C( I, J ) = BETA*C( I, J )
  214.    30             CONTINUE
  215.                   C( J, J ) = BETA*REAL( C( J, J ) )
  216.    40          CONTINUE
  217.             END IF
  218.          ELSE
  219.             IF( BETA.EQ.REAL( ZERO ) )THEN
  220.                DO 60, J = 1, N
  221.                   DO 50, I = J, N
  222.                      C( I, J ) = ZERO
  223.    50             CONTINUE
  224.    60          CONTINUE
  225.             ELSE
  226.                DO 80, J = 1, N
  227.                   C( J, J ) = BETA*REAL( C( J, J ) )
  228.                   DO 70, I = J + 1, N
  229.                      C( I, J ) = BETA*C( I, J )
  230.    70             CONTINUE
  231.    80          CONTINUE
  232.             END IF
  233.          END IF
  234.          RETURN
  235.       END IF
  236. *
  237. *     Start the operations.
  238. *
  239.       IF( LSAME( TRANS, 'N' ) )THEN
  240. *
  241. *        Form  C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) +
  242. *                   C.
  243. *
  244.          IF( UPPER )THEN
  245.             DO 130, J = 1, N
  246.                IF( BETA.EQ.REAL( ZERO ) )THEN
  247.                   DO 90, I = 1, J
  248.                      C( I, J ) = ZERO
  249.    90             CONTINUE
  250.                ELSE IF( BETA.NE.ONE )THEN
  251.                   DO 100, I = 1, J - 1
  252.                      C( I, J ) = BETA*C( I, J )
  253.   100             CONTINUE
  254.                   C( J, J ) = BETA*REAL( C( J, J ) )
  255.                END IF
  256.                DO 120, L = 1, K
  257.                   IF( ( A( J, L ).NE.ZERO ).OR.
  258.      $                ( B( J, L ).NE.ZERO )     )THEN
  259.                      TEMP1 = ALPHA*CONJG( B( J, L ) )
  260.                      TEMP2 = CONJG( ALPHA*A( J, L ) )
  261.                      DO 110, I = 1, J - 1
  262.                         C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
  263.      $                                          B( I, L )*TEMP2
  264.   110                CONTINUE
  265.                      C( J, J ) = REAL( C( J, J ) )         +
  266.      $                           REAL( A( J, L )*TEMP1 +
  267.      $                                 B( J, L )*TEMP2   )
  268.                   END IF
  269.   120          CONTINUE
  270.   130       CONTINUE
  271.          ELSE
  272.             DO 180, J = 1, N
  273.                IF( BETA.EQ.REAL( ZERO ) )THEN
  274.                   DO 140, I = J, N
  275.                      C( I, J ) = ZERO
  276.   140             CONTINUE
  277.                ELSE IF( BETA.NE.ONE )THEN
  278.                   DO 150, I = J + 1, N
  279.                      C( I, J ) = BETA*C( I, J )
  280.   150             CONTINUE
  281.                   C( J, J ) = BETA*REAL( C( J, J ) )
  282.                END IF
  283.                DO 170, L = 1, K
  284.                   IF( ( A( J, L ).NE.ZERO ).OR.
  285.      $                ( B( J, L ).NE.ZERO )     )THEN
  286.                      TEMP1 = ALPHA*CONJG( B( J, L ) )
  287.                      TEMP2 = CONJG( ALPHA*A( J, L ) )
  288.                      DO 160, I = J + 1, N
  289.                         C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
  290.      $                                          B( I, L )*TEMP2
  291.   160                CONTINUE
  292.                      C( J, J ) = REAL( C( J, J ) )         +
  293.      $                           REAL( A( J, L )*TEMP1 +
  294.      $                                 B( J, L )*TEMP2   )
  295.                   END IF
  296.   170          CONTINUE
  297.   180       CONTINUE
  298.          END IF
  299.       ELSE
  300. *
  301. *        Form  C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A +
  302. *                   C.
  303. *
  304.          IF( UPPER )THEN
  305.             DO 210, J = 1, N
  306.                DO 200, I = 1, J
  307.                   TEMP1 = ZERO
  308.                   TEMP2 = ZERO
  309.                   DO 190, L = 1, K
  310.                      TEMP1 = TEMP1 + CONJG( A( L, I ) )*B( L, J )
  311.                      TEMP2 = TEMP2 + CONJG( B( L, I ) )*A( L, J )
  312.   190             CONTINUE
  313.                   IF( I.EQ.J )THEN
  314.                      IF( BETA.EQ.REAL( ZERO ) )THEN
  315.                         C( J, J ) = REAL(        ALPHA  *TEMP1 +
  316.      $                                    CONJG( ALPHA )*TEMP2   )
  317.                      ELSE
  318.                         C( J, J ) = BETA*REAL( C( J, J ) )         +
  319.      $                              REAL(        ALPHA  *TEMP1 +
  320.      $                                    CONJG( ALPHA )*TEMP2   )
  321.                      END IF
  322.                   ELSE
  323.                      IF( BETA.EQ.REAL( ZERO ) )THEN
  324.                         C( I, J ) = ALPHA*TEMP1 + CONJG( ALPHA )*TEMP2
  325.                      ELSE
  326.                         C( I, J ) = BETA *C( I, J ) +
  327.      $                              ALPHA*TEMP1 + CONJG( ALPHA )*TEMP2
  328.                      END IF
  329.                   END IF
  330.   200          CONTINUE
  331.   210       CONTINUE
  332.          ELSE
  333.             DO 240, J = 1, N
  334.                DO 230, I = J, N
  335.                   TEMP1 = ZERO
  336.                   TEMP2 = ZERO
  337.                   DO 220, L = 1, K
  338.                      TEMP1 = TEMP1 + CONJG( A( L, I ) )*B( L, J )
  339.                      TEMP2 = TEMP2 + CONJG( B( L, I ) )*A( L, J )
  340.   220             CONTINUE
  341.                   IF( I.EQ.J )THEN
  342.                      IF( BETA.EQ.REAL( ZERO ) )THEN
  343.                         C( J, J ) = REAL(        ALPHA  *TEMP1 +
  344.      $                                    CONJG( ALPHA )*TEMP2   )
  345.                      ELSE
  346.                         C( J, J ) = BETA*REAL( C( J, J ) )         +
  347.      $                              REAL(        ALPHA  *TEMP1 +
  348.      $                                    CONJG( ALPHA )*TEMP2   )
  349.                      END IF
  350.                   ELSE
  351.                      IF( BETA.EQ.REAL( ZERO ) )THEN
  352.                         C( I, J ) = ALPHA*TEMP1 + CONJG( ALPHA )*TEMP2
  353.                      ELSE
  354.                         C( I, J ) = BETA *C( I, J ) +
  355.      $                              ALPHA*TEMP1 + CONJG( ALPHA )*TEMP2
  356.                      END IF
  357.                   END IF
  358.   230          CONTINUE
  359.   240       CONTINUE
  360.          END IF
  361.       END IF
  362. *
  363.       RETURN
  364. *
  365. *     End of CHER2K.
  366. *
  367.       END
  368.